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Avalanches of grain displacements can be generated by creating local voids within the interior of a 
granular material at rest in a bin. Modeling such a two-dimensional granular system by a collection 
of mono-disperse discs, the system on repeated perturbations, shows all signatures of Self-Organized 
Criticality. During the propagation of avalanches the competition among grains creates arches and 
in the critical state a distribution of arches of different sizes is obtained. Using a cellular automata 
model we demonstrate that the existence of arches determines the universal behaviour of the model 
system. 

PACS numbers: 05.70.Jk, 64.60 Lx, 74.80.Bj, 46.10,+z 



o 

CO 



I 

c 

o 
o 



> 
o 

CO 
O 
OO 

o 

OO 



T3 

c 

O 
o 



The search for self organized criticality (SOC) [1] in 
granular systems in general and in sandpiles in particular 
has been a subject of active research over the last decade. 
How such a system reacts in the form of fast cascades of 
grain displacements, called 'avalanches', in response to a 
slow external drive at the microscopic level is the cru- 
cial question of study. It has been suggested that start- 
ing from an arbitrary initial condition, a non-equilibrium 
critical state, characterised by scale free avalanches in 
both space and time, should be expected after a long 
time [1,2]. Other naturally occuring physical phenomena 
like forest fires [3], river networks [4], earthquakes [5] etc. 
have also been proposed as examples of systems showing 
SOC. 

It was observed that a steady shaped sandpile grown 
on a fixed base fulfills all the requirements of SOC. In this 
state, dropping even a few grains creates rapidly moving 
avalanches of sand sliding along the surface of the pile. 
It was expected that the avalanches should be power law 
distributed in their spatial as well as temporal extents 
and therefore a sandpile should be considered as a simple 
example showing SOC [1,2]. 

Experimental observations, however, show partial sup- 
port to this idea. Sand is allowed to fall from a slowly ro- 
tating semi-cylindrical drum through the space between 
the plates of a vertical parallel plate capacitor. The 
Fourier spectrum of the time series data of fluctuating 
capacitance showed a peak contrary to the expectation 
of a power law [6] . Similarly sand dropping from the edge 
of a sandpile on a fixed base was directly measured, and 
was seen to have a power law distribution only for small 
systems but not for large systems [7] . It was argued that 
due to the approximately spherical shape of the grains 
used in these experiments the effect of inertia cannot be 
neglected and this was the reason for absence of scaling 
behaviour. This is verified in an experiment using rice 
grains, which are highly anisotropic, and criticality was 
observed [8]. 

A number of theoretical models, generally known as 



'sandpile' models have been studied. The models are 
based on stochastically driven cellular automata which 
evolve under a nonlinear, diffusive mechanism leading to 
a nonequilibrium critical state. The most widely studied 
is the Abelian sandpile model (ASM) where the stable 
configuration does not depend on the sequence of sand 
grain additions [9] . Other variants of the sandpile model 
include situations where the stability of a sand column 
depends on the local slope or the Laplacian of the height 
profile [10]. A two-state sandpile model with stochastic 
evolution rules is also studied [11]. 

In all the studies discussed above the avalanches prop- 
agate on the surfaces of the sandpiles. However, there 
exists the possibility of creating avalanches in the inte- 
rior of a granular material. In a granular material kept in 
a bin at rest, different grains support one another by mu- 
tually acting balanced forces. Now if a grain is removed, 
the grains which were supported by it become unstable 
and tend to move. Eventually the grains in the further 
neighbourhood also loose their stability. As a result an 
avalanche of grain displacements takes place, which stops 
when no more grains remain unstable. The basic physical 
behavior here is thus quite different from the avalanches 
on the surface of the pile because of the constraints to 
particle motion in the dense particle beds. 

A two-dimensional semi-lattice model was studied for 
this problem [12]. Non-overlaping unit square boxes 
model the grains, whose horizontal coordinates can vary 
continuously where as the vertical co-ordinates are dis- 
cretized. A grain can only fall vertically if insufficiently 
supported and sufficient space below is available. The 
system is disturbed by removing grains at the bottom 
and thus creating avalanches of grain movements. 

During the propagation of avalanches inside a granular 
medium grains compete locally with one another to oc- 
cupy the same vacant space. The high packing densities 
of the particle beds prevent a single particle from occupy- 
ing the available void space, and consequently particles 
get locked to form 'arches' [13]. A stable arch is a chain 
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of grains where the weight of each grain is balanced by 
the reaction forces from two neighbouring grains in the 
chain. Arches can form only when a grain is allowed to 
roll over other grains. Since the rolling motion of the 
grains was absent therefore the arches were not formed 
in the granular patterns [12]. 

We here study a more realistic model of this problem 
in two dimension, where both fall and roll motions of 
grains are allowed. The granular system is represented 
by N hard monodisperse discs of radii R. No two grains 
are allowed to overlap but can touch each other and one 
can roll over the other if possible without slipping. A 
rectangular area of size L x L on the x—y plane represents 
our two dimensional bin. Periodic boundary condition is 
used along the x direction and the gravity acts along the 
—y direction. The bottom of the bin at y = is highly 
sticky and any grain which comes in contact gets stuck 
there and does not move further. 

The dynamical evolution of the system is studied by a 
'pseudo-dynamics' [14]. Unlike the method of molecular 
dynamics we donot solve here the classical equations of 
motion for the grain system. Only the direction of gravity 
and the local geometrical constraints due to the presence 
of other grains govern the movement of a grain. For 
justification of using the pseudo-dynamics we argue that 
due to the high compactness of the system a grain never 
gets sufficient time to accelerate much. Therefore in our 
model, a grain can have only two types of movements in 
unit time: The vertical fall up to a maximum distance 
6 and the roll up to a maximum angle 9 = S/2R over 
another grain in contact. 

Movement of a grain needs the information on the po- 
sitions of other grains in the neighbourhood which it may 
possibly interact. An efficient way to search this is to dig- 
itize the bin into a square grid and associate the serial 
number of a grain to the primitive cell of the grid con- 
taining its centre. Choise of R = l/y/2+ ensures that a 
cell corresponds to at most one grain. With this choice 
it is sufficient to search only within the 24 neighbouring 
cells for possible contact grains. The weight of a grain n 
is supported by two other grains. If til and ur are the 
serial numbers of the left (L) and right (R) supporting 
grains, then the grain n is updated as: (i) If til = ur = 0, 
it falls, (ii) if nj, ^ but ur = 0, it rolls over ul, (hi) if 
Ul = but riR ^ 0, it rolls over ur, (iv) if ul i= and 
ur it is stable. 

When the grain n with the centre at (x n , y n ) is allowed 
to fall, it may come in contact with a grain r at (x r , y r ) 
within the distance S. The new coordinates are: 



coordinates for the centre of n where it touches both r 
and t are: 



Vn 



Vr 



+ ^/AR 2 - (x n - x r ) 2 if y n -y' n <5 



otherwise, y' n = y n - 5. 

Similarly the grain n, while rolling over the grain r, may 
come in contact with another grain t at (xt, yt)- The new 
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Here d r t is the distance between the centres of the grains 
r and t and the g = +1 and -1 for left and right rolls. 
To reach the new position if the grain n rolls an angle 
9 m < 9 then it is accepted, otherwise it rolls up to an 
angle 9. 

The initial grain pattern is generated by the ballistic 
deposition method with restructuring [15,16]. Grains are 
released sequentially one after another along randomly 
positioned vertical trajectories and are allowed to fall till 
they touch the growing pile and then roll down along the 
paths of steepest descent to their local stable positions. 
We notice that in the initial pattern no big arch exists 
since a grain while rolling down along the surface does not 
need to compete with any other grain. Using a system 
of L=80, we compute the packing fraction p = 0.822 ± 
0.005 consistent with more precise estimate of 0.8180 ± 
0.0002 in [16]. 




FIG. 1. Picture of a typical avalanche. The vacant cir- 
cles denote the undisturbed grains, filled circles denote the 
displaced grains, the opaque circle at the bottom denotes the 
position of the grain which was removed and the shaded cir- 
cle denote the position of the grain where the removed grain 
was replaced. In a system of size L=30 and number of grains 
iV=340, 97 grains took part in the avalanche. 

The system is repeatedly perturbed by removing grains 
randomly at the bottom one after the other. Every time 
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a grain is removed, an avalanche is followed and after 
the avalanche is over the removed grain is placed back 
at a random position on the top surface again using the 
same ballistic method (Fig. 1). We first observe that 
p of the stable configuration, averaged over many initial 
random patterns generated with the same algorithm, de- 
creases with the number of avalanches created and finally 
reaches a steady value of 0.748 ± 0.005. A similar study 
but with different initial configurations with a different 
value of average initial p, shows the final steady state 
packing fraction reaches the same value, which implies 
that the final stable state is independent of the initial 
state. The total number of grain displacements is called 
the avalanche size s and the duration of the avalanche is 
the life time t. Power law distributions are observed for 
s and t: D(s) ~ s" 7 " 3 and D{t) ~ i~ Tt with r s w 1.7 (Fig. 
2) and r t w 2.1. 
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FIG. 2. Plot of the avalanche size distribution D(s) vs. 
s for a system of size N = 10000 and over 70000 avalanches. 
The straight part fits with a slope t s = 1.7. 

The observed exponents in this model are significantly 
larger than the corresponding values of t s w 1.34 and 
tl ~ 1-47 in [12]. We argue that the reason for this 
difference could be the existence of arches in our model. 
Due to the arches, the propagation of avalanches get ar- 
rested more frequently than in [12] and therefore smaller 
avalanches are more probable which enhances the values 
of the critical exponents in our model. 

To demonstrate more explicitly that the above reason- 
ing may be correct, we study a cellular automata model 
for the granular systems. A square lattice of size L with 
periodic boundary condition along the x axis represents 
the bin. The gravity acts in the — y direction. Positions 
of the grains are limited only to the lattice sites: a site 
can be either occupied by a grain or remain vacant. Ini- 
tially the bin is filled up again using the ballistic method. 

A single movement of a grain at C(i,j) involves 



the neighbouring seven sites: LU(i — 1, j + l),L{i — 
l,j),LD(i-l,j-l),D(i,j-l),RD(i+l,j-l),R(i+l,j) 
and RU(i + 1, j + l).In the fall move the grain comes 
down one level to the vacant site at D and in roll move 
the grain goes either to the vacant site at LD or at RD 
(Fig. 3). An arch is formed when a grain at C is con- 
sidered stable if both of its two diagonally opposite sites 
either at LD and RU or, at LU and RD are occupied. As 
a result, on the square lattice, the only possible shape of 
an arch consists of two sides of a triangle. Depending on 
whether we allow the arch formations or not, we define 
the following two models. 
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FIG. 3. The possible fall and roll moves in the cellular 
automata model of the granular system on the square lattice. 
Filled circle denotes the position of a grain, unfilled circle 
denotes a vacant site. The grain moves to the vacant position 
irrespective of the occupation of the sites with square boxes. 
Shaded rectangle denotes a pair of sites, in which at least one 
is occupied. 

In model A we allow arch formations. The grain at 
C falls only if any of the following three conditions is 
satisfied: (i) LD, L and LU are vacant (ii) RD, R and 
RU are vacant (iii) LD and RD are vacant. In all other 
situations the grain does not fall. Notice that in condi- 
tions (i) and (ii) we are allowing the formation of arches. 
The grain at C rolls only if the site D is occupied. This 
is done in any of the three following ways: (i) If both 
LD, L are vacant but either of RD, R is occupied then 
the grain rolls to LD. (ii) if both RD, R are vacant but 
either of LD, L is occupied then the grain rolls to RD. 
(iii) If all four sites at LD, L, RD, R are vacant then the 
grain rolls either to LD or to RD with probability 1/2. 
A steady state pattern is shown in Fig. 4. 

In model B we do not allow arch formations. The first 
two conditions for fall of model A are modified as: (i) LD 
and L are vacant (ii) RD and R are vacant. All other 
conditions of fall as well as roll remain same as in model 
A. 
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Initial granular patterns are generated using the same 
random ballistic deposition method with restructuring 
[15,16] and patterns are same for both models. As before 
the avalanches are created by taking out one grain at 
a time at the bottom, allowing the system to relax and 
replacing it randomly on the surface after the avalanche 
is over. Here also we see that the average density of sites, 
starting from an initial value of 0.907 ± 0.005, decreases 
to the final stable value of 0.590±0.005 in model A and to 
0.618±0.005 in model B. The avalnche size s and life time 
t follow power law distributions: D(s) <~ s~ r s , D(t) ~ 
t~ Tt and similarly for the model B. Different exponents 
are obtained for the two models: w 1.48 and s=s 
1.99 where as rf w 1.34 and rf w 1.50. We explain 
that absence of arches makes the exponent values for the 
model B close to that of [12] but their presence enhances 
the values in model A. 



FIG. 4. Figure 4: A stable configuration of the granular 
system in the cellular automata model A after a large number 
of avalanches are created. Triangular arches are noticed. 

To summarise, avalanches of grain displacements can 
be created in the interior of a granular material at rest 
by locally disturbing the system. In a numerical study 
we provide indications that on repeated creations of such 
avalanches, the granular system reaches a critical state 
characterised by long range correlations. The presence 
of arches play an important role in determining the crit- 
ical behaviour. Since in an avalanche grain motions are 
highly constrained, the effect of inertia may not be very 
significant. Therefore, we conjecture that in real experi- 
ment of internal avalnches it should be possible to observe 
SOC even with spherical grains unlike the case of surface 



avalanches, where anisotropic grains were necessary to 
observe criticality. 

The funding provided by the Indo-French Centre for 
the Promotion of Advanced Research (IFCPAR) is grate- 
fully acknowledged. 

[1 ] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. 
Lett. 59. 381 (1987). 

[2 ] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. 
A 38, 364 (1988); P. Bak, How Nature Works: The 
Science of Self-Organized Criticality, (Copernicus, 
New York, 1996). 

[3 ] P. Bak and K. Chen, Physica D 38, 5 (1989), H- 
M Broker and P. Grassberger, Phys. Rev. E 56, 
R4918 (1997). 

[4 ] A. Rinaldo, I. Rodriguez-Iturbe, R. Rigon, E. 
Ijjasz-Vasquez and R. L. Bras, Phys. Rev. Lett. 
70, 822 (1993); S. S. Manna and B. Subramanian, 
Phys. Rev. Lett. 76 (1996) 3460. 

[5 ] J. M. Carlson and J. S. Langer, Phys. Rev. Lett. 
62, 2632 (1989); Z. Olami, H. J. S. Feder and K. 
Christensen, Phys. Rev. Lett. 68, 1244 (1992). 

[6 ] H. M. Jaeger, C-h Liu and S. R. Nagel, Phys. Rev. 
Lett. 62, 40 (1989). 

[7 ] G. A. Held, D. H. Solina II, D. T. Kcane, W. J. 
Haag, P. M. Horn and G. Grinstein, Phys. Rev. 
Lett. 65, 1120 (1990). 

[8 ] V. Frette, K. Christensen, A. Malthe-Sorensscn, J. 
Feder, T. Jossang and P. Meakin, Nature (London) 
379, 49 (1996). 

[9 ] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990). 

[10 ] L. P. Kadanoff, S. R. Nagel, L. Wu and S. Zhou, 
Phys. Rev. A. 39, 6524 (1989); S. S. Manna, Phys- 
ica A 179, 249 (1991). 

[11 ] S. S. Manna, J. Phys. A 24, L363 (1992). 

[12 ] R. E. Snyder and R. C. Ball, Phys. Rev. E 49, 
104 (1994). 

[13 ] D. E. Wolf, in Computational Physics ed. by K. 
H. Hoffmann and M. Schrcibcr, Springer 1996. 

[14 ] S. S. Manna and D. V. Khakhar, in Nonlin- 
ear Phenomena in Material Science III, G. Anan- 
thakrishna, L. P. Kubin and G. Martin (Eds.) 
(Transtech, Switzerland, 1995). 

[15 ] W. M. Visser and M. Bolsterli, Nature 239, 504 
(1972). 

[16 ] R. Jullien, P. Meakin and A. Pavlovitch, in Disor- 
der and Granular Media, D. Bideu and A. Hansen 
(Eds.) (Elsevier 1993). 




4 



